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SUMMARY 


Numerical  FORTRAN  algorithms  are  presented  for  cxth  quantiles  cf  the 
distribution  of  an  i.i.d.  random  variable  distributed  according  to  a  hyper-Gamma,  or  a 
generalized  Gumbel,  or  a  lognormal  probability  density  function.  They  are  based  on 
highly  efficient  routines  for  the  evaluation  of  the  incomplete  Gamma  and  error 
functions.  A  program  diskette  will  be  made  available  upon  request. 
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INTRODUCTION 


To  obtain  information  about  a  distribution  it  is  frequently  desirable  to  knew  its 
quantiles  for  certain  levels  o<  .  In  general  terms,  let  X  be  an  i.i.d.  random  variable 
with  continuous  probability  density  function  (pdf)  f(x;P)  where  P  denotes  a 
parameter  vector  and  -«  <  x  <  ». 

Let 


x 

F(x)  =  Pr  (x  s;  x)  =  J  f  (t ; P)  dt 

-00 


be  the  associated  cumulative  distribution  function  (cdf).  The  «th  quantile  (ren- 
exceedance  probability  at  level  cx  )  of  the  distribution  of  X  is  that  value  x,*  for 
which 


Ffcoc) 


CX  . 


<") 


It  is  also  referred  to  as  the  (lOOcx)th  percentile  of  the  distribution  of  X  .  Since  0  < 
F(x)  <  I  and  since  F(x)  is  strictly  monotonically  increasing  it  is  evident  that  <x  must 
satisfy  the  inequality  0  <  <x  <  1  and  that  equation  (*)  has  exactly  one  solution. 

Equation  (*)  may  also  be  expressed  in  the  form  x*  =  G(o<)  where  G  is  the  inverse  of  F. 

Unfortunately,  there  are  only  a  few  distributions  for  which  the  cdf  can  be 
inverted  in  closed  form,  a  particular  one  for  which  this  is  possible  is  the  exponential 
distribution  with  pdf  f(x;A,b)  =  b"1  exp  -£  ,  £  =  (x-A)b-'.  In  general,  therefore,  it  is 
necessary  to  solve  the  quantile  equation  (")  numerically  for  given  levels  cx  once  the 
distribution  pdf  f(x;P)  has  been  specified. 

This  report  presents  efficient  algorithms  for  the  numerical  solution  of  the 
quantile  equation  (*)  for  three  specific  classes  of  probability  distributions: 

(i)  Hyper-Gamma: 


f(x;P) 


f,  -  £"p  exp  -£^  ,  £  =  (x-X)P_1  ,  x  >  A 
bT(a) 

k  0  ,  x  <  A  , 


0-0 


a  =  ( J -p)3 ” 1  ,  P  =  (A,  b,  3,  p),  where  A  is  the  real  shift  parameter,  b  >  0  the  scale 
parameter,  3  >  0  the  terminal  shape  parameter,  and  p  <  I  the  initial  shape  parameter. 
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(ii)  Generalized  Gumbel  (Extreme  Value  Type  I  for  Maximum  Elements): 


f(x;P)  =  .  ^  exp  -  £  f  e  ^  -  s)  .  £  =  (x-X)b-1  ,  -oo  <  x  <  «  ,  (1.2) 

t>r(3)  v  J 

P  =  (X,  b,  £),  where  X  is  the  real  shift,  b  >  0  is  scale,  and  3  >  0  is  shape.  For  £  =  1 
the  classical  Gumbel  pdf  results. 


P  =  (X.a.p),  X  is  real  shift,  c  >  0  is  shape,  is  real  scale.  (The  scaling  property  c 
jj  becomes  apparent  if  one  sets  jj  =  log  b,  b  >  0.) 

For  the  three  distribution  classes  defined  by  the  pdf’s  (1.1),  (1.2),  and  (1.3), 
parameter  estimation  algorithms  for  given  data  have  recently  been  made  available  [l], 
[2,3],  [4].  Their  parameter  vector  outputs  could  directly  be  used  as  inputs  for  the 
quantile  algorithms  offered  in  this  report  in  order  to  achieve  a  fully  automated 
distribution  assessment  package  if  it  were  augmented  by  a  goodness-of-f it  evaluation 
algorithm.  However,  goodness-of-f  it  considerations  have  been  beyond  the  scope  of  the 
investigations  that  resulted  in  the  parameter  estimation  methods,  and  they  are  outside 
the  scope  of  the  present  investigation.  Furthermore,  different  analysts  may  have 
different  objectives  depending  on  the  ultimate  use  of  the  numerical  results. 

Therefore,  the  three  quantile  algorithms  are  presented  in  stand-alone  versions. 
To  facilitate  easy  hook-up  to  the  existing  parameter  estimation  algorithms  they  have 
been  packaged  individually.  A  program  diskette  will  be  made  available  upon  request. 
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11. 


THE  QUANTILE  EQUATIONS 


Throughout  this  section  let  y  =  x  -  X,  and  r\  =  yb"1. 

(i)  The  cdf  for  the  hyper -Gamma  distribution  class  (1.1)  takes  the  form 

F (y)  =  •— rr  *(aV)  .  i\  >0  . 
i  (a)  v  J 

with  a  =  (l-p)fr1  ,  r(x)  the  (complete)  Gamma  function,  and  Z(v,x)  the  incomplete 
Gamma  function  with  upper  integration  limit  x  . 

With  =  u  >  0  the  general  quantile  equation  (*)  reduces  to  2f(a,u)  -  c<r(a)  =  C. 
For  computational  efficiency  it  is  useful  to  transform  this  equation  into  the  mere 
convenient  form 

9(u)  =  tf*(a,u)  -  <x  exp  -  (a  log  u)  =  C  ,  «  e  (0, ')  ,  (2.0 

where  Z*  (a,u)  is  the  reduced  incomplete  Gamma  function,  Z*  (a,u)  =  2f(a,u)/[uar(a)]. 

If  Uo<  is  the  (unique)  solution  of  (2.1),  the  corresponding  x  value  is  given  by 

x^  =  X  -  bexp  ^3  1  log  .  (2.1a) 

(ii)  The  cdf  for  the  generalized  Gumbel  distribution  class  (1.2)  is 

F<g)  =  rk 

with  3  >  0,  T(v,x)  the  complementary  incomplete  Gamma  function,  r(v,x)  =  r(v)  -  tf(v,x). 

With  3  exp  -T[  =  u  >  0  equation  (*)  now  takes  the  form  r(3,u)  -  o<r(S)  =  0  which 
can  be  changed  into  the  more  convenient  form 

<P(u)  =  tf“(3.u)  -  (l-<x)  exp  -  (3  log  u)  =  C  ,  <x  e  (0.0  .  (2.2) 

If  Uo<  is  the  (unique)  solution,  one  obtains  the  corresponding  x  value  as 


Xo<  =  X  -  b  log  (3  u*'1) 


(2.2a) 
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(lii)  Under  the  transformation  (a  J~2)  1  leg  r  =  u  e  ,  (b  =  exp  p),  the  cdf  for 
the  lognormal  distribution  class  (1.3)  can  be  written  as 

F (y)  =  1(1  ‘erf  u) 

with 

u 

erf  u  =  J  exp  -  t^dt  . 
o 

The  quantile  equation  to  be  solved  then  is 

«P(u)  =  erf  u  -  (2cx  -  ’)  =  0  .  <x  e  (0,1)  .  (2.3) 


If  Uo<  is  the  (unique)  solution,  the  corresponding  x  value  is  given  by 

x^  =  X  -  exp  (p*a  .  (2.3a) 
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DESCRIPTION  OF  THE  ALGORITHMS 


In  each  of  the  three  cases  discussed  here  the  user  has  tc  provide  as  mput  a 
complete  set  of  admissible  distribution  parameter  values.  The  user  t her  has  the  option 

either  to  set  o<  at  any  desired  level  in  the  interval  (0,1)  or  to  ask  for  quantiles  for  a 

pre-established  set  of  fairly  standard  c<  values.  (The  program  printout  should  be 
consulted  for  details.) 

(i)  The  solution  algorithms  for  equation  (2.1)  is  started  by  evaluation  of  <?(l). 
If  ip(i)  =  c,  then  uc  =  1.  and  x0  =  X  -  b  according  tc  (2.1a).  If  <P(1)  <  C,  the  root  uc  of 
(2.1)  is  greater  than  unity.  In  this  case  the  function  <P(u)  is  evaluated  at  the  points  of 
the  sequence  uv  =  1  *  1C"1  (1.6'3)v_1  (v=1.  2,  ...)  until  the  root  u*  is  bracketed  by 

change  in  sign  of  <p  .  If  f(1)  >  0,  the  root  uw  t  (0,1).  New  the  search  sequence  uv  = 

1  -  2V-1/(2V"'  *  99)  (v=l,  2.  ...)  is  used  tc  achieve  root  bracketing. 

Cnee  Uc*  has  been  bracketed,  Brent's  method  [5;  Chap.  7.3]  performs  final 
calculation  of  u,*  .  The  algorithm  returns  by  (2.1a). 

The  function  £*(a,u)  is  evaluated  by  means  of  a  routine  provided  in  [6;  Chap. 

45.3]. 


(n)  Because  of  the  structural  similarity  between  equations  (2.1)  and  (2.2)  the 
solution  strategy  employed  relative  tc  the  first  carries  ever  directly  tc  the  ether.  The 
only  difference  is  that  new  fO)  >  0  implies  that  u<x  >  1,  <P(1)  <  0  implies  that  e 
(0,1). 


(lii)  The  solution  algorithm  of  equation  (2.3)  is  initiated  by  evaluation  of  <f(G). 

If  cx  =  1/2,  then  ui/2  =  0,  and  xi/2  =  X  *  exp  jj  according  to  (2.3a).  If  <f>(0)  <  C,  (1/2 
<  cx  <  1).  the  root  u<*  >  0.  The  search  sequence  uv  =  10"1  (1.618)V_1  (v=1,2,  ...)  is 
used  for  root  bracketing.  If  f(0)  >  0  (0  <  cx  <  1/2),  Uc*  <  C.  The  search  sequence 
uv  =  -  1 0” 1  (1.61  Q) v“ 1  (v=i,  2,  ...)  is  used  now.  Brent's  method  is  employed  again. 

Error  function  evaluation  is  dene  by  a  routine  given  in  [6;  Chap.  40.8]. 

Warning:  Because  of  properties  of  the  functions  %*  and  erf  it  may  happen 

that  either  of  the  three  quantile  algorithms  is  unable  to  produce  a  root  of  the 
corresponding  equation  <p(u)  =  0.  This  is  true  if  the  parameters  a  and  0  in  equations 
(2.1)  and  (2.2),  respectively,  are  small  or  large.  Relative  to  equation  (2.3)  failure  may 
occur  if  o<  or  1-cx  is  small. 


\ 
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IV,  EXAMPLES 

Example  1. 

Three-parameter  hyper-Gamma  estimation  (with  X  assumed  tc  be  zero)  has  Deer 
performed  in  (I;  Sec.  7,  Ex.  2,  Table  2.1]  for  (simulated)  observations  of  interarrival 
times  (in  minutes)  at  a  drive-up  banking  facility  over  a  90-minute  period  [7;  Chap.  5.3  1, 
Ex.  5.1,  p.  1  84],  The  smallest  and  largest  observations  are  Xmin  =  0.1,  xmax  =  1.96. 
Cuantiles  ter  various  cx  levels  for  the  moment  (n)  and  maximum-likeliheed  (ML) 
parameters  are  shewn  in  Table  1  of  Sec.  V.  The  x^  values  are  in  minutes. 

LW-B-L.Q-2. 

Four-parameter  hyper-Gamma  estimation  has  beer  performed  on  another  data  set 
[1:  Sec.  7,  Ex.  3,  Table  3.2]  which  represents  observations  of  daily  mean  temperatures 
(in  °F)  at  Shemya,  Alaska,  during  the  winters  of  1960-1979.  with  xmin  =  18,  xmax  = 

40.  The  results  for  M  and  ML  parameters  are  displayed  ir  Table  2.  The  x^  values  are 
in  °F. 


L\arr.r-le  2  • 

Three-parameter  generalized  G umbel  estimation  results  for  observations  of 
annual  24-hour  maximum  rainfalls  (in  points,  xmin  =  1-54,  xmax  =  1  105)  at  Sidney, 
Australia,  over  the  period  1859-1945,  have  been  provided  in  [3;  Sec.  7,  Ex.  1,  Table  2], 
The  original  data  are  from  [8],  Table  3  shows  the  M  and  ML  quantiles.  The 
values  are  in  points. 

Example  4- 

Another  three-parameter  generalized  Gumbel  estimation  has  been  done  in  [3;  Sec. 
7,  Ex.  2,  Table  5]  for  observations  of  annual  peak  gust  winds  [in  kts]  at  Argentia  Naval 
Air  Station,  Newfoundland,  over  the  period  1941-1963  (original  data  from  [9;  Chap. 
4.12],  xmin  =  62,  xmax  =  91).  Results  are  shown  in  Table  4;  xw  is  in  kts. 


Eaacn&lg  5- 

Two-parameter  lognormal  estimation  (X  assumed  to  be  zero)  has  been  performed 
in  [4;  Sec.  VIII,  Ex.  4,  Table  4.3]  on  observations  of  rainfall  totals  (in  inch,  xmin  =  2.5, 
xmax  =  18.5)  for  sets  of  four  consecutive  months  at  Camden  Square,  London,  over  the 
period  1870-1943  [10;  Chap.  7.5,  Ex.  7.51  1].  Table  3  shows  the  results;  x^  in  inch. 
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The  report  [4;  Sec.  VIII,  Ex.  5,  Table  5.31  contains  another  two-parameter  (X=C) 
lognormal  estimation  for  data  given  in  [I  1;  Ex.  1C,  p.  47).  They  represent  weekly 
precipitation  sums  observations  (in  10~2  inch,  xmin  =  25,  xmax  =  825)  at  Kwajalem, 
Marshall  Islands,  during  the  summers  1949-1958.  Quantile  results  are  displayed  in 
Table  6;  x,*  is  in  1C-2  inch. 


7, '(8  Blank) 
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TABLE  1.  M  and  ML  Quantiles,  Drive-Up  Bank 


DISTRIB.TYPE:  Pour  Parameter  Hyper-Gamma 


LAMBDA  - 
BETA 
BVAL 
PVAL 


Mom.  Estim. 


0.000000 

1.625230 

0.809469 

0.265411 


ML.  Estim. 


0.000000 

0.967295 

0.354808 

-0.069105 


Cum.Pct.Lvl.  Quantiles 

0.0102 

0.000002 

0.000067 

0.1002 

0.000057 

0.000581 

1.0002 

0.001300 

0.005038 

2.0002 

0.003339 

0.009703 

5.0002 

0.011628 

0.023324 

10.0002 

0.029920 

0.046070 

25.0002 

0.105559 

0.119964 

50.0002 

0.288488 

0.280036 

75.0002 

0.583231 

0.550817 

90.0002 

0.927675 

0.908438 

95.0002 

1.161985 

1.179608 

98.0002 

1.446185 

1.539171 

99.0002 

1.646047 

1.812012 

99.9002 

2.243029 

2.723128 

99.9902 

2.768327 

3.639656 

TABLE  2. 

M  and  ML  QuaQtiles 

,  Shemya 

DISTRIB.TYPE:  Four  Parameter  Hyper-Gamma 

Mom.  Estim. 

ML.  Estim. 

LAMBDA  « 

12.750255 

13.517533 

BETA 

7.575647 

7.952291 

BVAL 

22.152744 

21.772694 

PVAL 

-3.979800 

-3.593978 

Cum.Pct.Lvl.  Quantiles 

0.0102 

16.163086 

16.377214 

0.1002 

18.169350 

18.238099 

1.0002 

21.355516 

21.310090 

2.0002 

22.641807 

22.579652 

5.0002 

24.646488 

24.583222 

10.0002 

26.441695 

26.396393 

25.0002 

29.317097 

29.317104 

50.0002 

32.164208 

32.196963 

75.0002 

34.593508 

34.620922 

90.0002 

36.462796 

36.458289 

95.0002 

37.461554 

37.429782 

98.0002 

38.495249 

38.428114 

99.0002 

39.137810 

39.045228 

99.9002 

40.772239 

40.603981 

99.9902 

41.966407 

41.733911 
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TABLE  3.  M  and  ML  Quantiles,  Sidney 
DISTRIB.TYPE:  Three  Parameter  Genl.Gumbel 


Mom.  Estim. 

ML.  Estim. 

LAMBDA  - 

348.336000 

335.783000 

BETA 

1.056140 

0.735367 

BVAL 

143.341000 

110.752000 

Cum.Pct .Lvl. 


Quantiles 


0.0102 

0.100Z 

1.0002 

2.000Z 

5.0002 

10.000Z 

25.0002 

50.0002 

75.0002 

90.0002 

95.0002 

98.0002 

99.0002 

99.9002 

99.9902 


35.486434 

76.136252 

133.353174 

156.320799 

193.824197 

230.716691 

301.590770 

397.861599 

519.170057 

656.844906 

755.103127 

881.978195 

976.912277 

1290.236179 

1602.838134 


66.041079 

100.207878 

149.118178 

169.072924 

202.093947 

235.141751 

300.316124 

392.554633 

514.573585 

659.252220 

765.194247 

903.887393 

1008.450385 

1355.340981 

1702.132745 
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TABLE  4.  M  and  ML  Quantiles,  Argentia 

DISTRIB.TYPE:  Three  Parameter  Genl.Gumbel 

Mom.  Estim.  ML.  Estim. 


LAMBDA  - 
BETA  - 
BVAL 


Cum.Pct.Lvl. 


0.010Z 

0.100Z 

l.OOOZ 

2.000Z 

5.000Z 

10.000Z 

25.000Z 

50.000Z 

75.000Z 

90.000Z 

95.000Z 

98.000Z 

99.000Z 

99.900Z 

99.990Z 


74.468400 

4.410460 

13.390300 


56.646226 

59.210458 

62.609585 

63.907460 

65.947261 

67.860895 

71.294876 

75.504766 

80.199135 

84.907464 

87.973928 

91.672466 

94.295208 

102.384949 

109.954114 


74.528300 

4.757460 

13.941600 


56.499598 

59.105693 

62.551503 

63.864472 

65.924875 

67.854326 

71.308055 

75.526940 

80.211754 

84.890403 

87.927217 

91.579797 

94.163588 

102.103430 

109.497965 


Quantiles 
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TABLE  5.  M  and  ML  Quantiles,  Camden  Square 


DISTRJB.TYPE:  Two-Parameter  Log' 


Mom.  Estim.  ML 


SIGMA  -  0.321706  0 

MU  -  2.085390  2 


Cua.Pct.Lvl.  Quantiles 

0.010Z 

2.432599 

2 

0.100Z 

2.977976 

2 

1.000Z 

3.807565 

3 

2.000Z 

4.156555 

3 

5.000Z 

4.740918 

4 

10.000Z 

5.328699 

5 

25.000Z 

6.477940 

6 

50.000Z 

8.047729 

8 

75.000Z 

9.997923 

10 

90.000Z 

12.154176 

12 

95.000Z 

13.661058 

14 

98.000Z 

15.581641 

16 

99.000Z 

17.009807 

17 

99.900Z 

21.748315 

23 

99.990Z 

26.624183 

28 

Normal 


Estim. 


341545 

081320 


250401 

789512 

621063 

974394 

570070 

173826 

365873 

015042 

091451 

416516 

056874 

163695 

740896 

029440 

546419 


TABLE  6.  M  and  ML  Quantiles,  Kwajalein 


DISTRIB.TYPE: 

Two-Parameter 

Log-Normal 

Mob.  Estim. 

ML.  Estim. 

SIGMA  » 

0.677532 

0.819763 

MU 

5.211890 

5.1A3320 

Cua.Pct.Lvl. 

Quantiles 

0.010Z: 

14.763196 

8.122255 

0.100Z: 

22.604720 

13.599893 

1.000Z: 

37.928958 

25.438559 

2.000Z: 

45.622882 

31.808471 

5.000Z: 

60.186557 

44.475116 

10.000Z: 

76.984015 

59.904489 

25.000Z: 

116.152224 

98.533664 

50.000Z: 

183.440433 

171.283487 

75.000Z: 

289.709411 

297.746288 

90.000Z: 

437.108827 

489.746816 

95.000Z: 

559.101469 

659.650503 

98.000Z: 

737.577083 

922.333953 

99.000Z: 

887.195286 

1153.289884 

99.900Z: 

1488.644532 

2157.225330 

99.990Z: 

2279.343398 

3612.054976 
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